#include "data.h"
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <omp.h>
#include <windows.h>

#define QUAD 9
#define MISSTAKE 1E-15

using namespace std;

void f( double r, double &force, double &potential )
{
	double ri = 1/r;
	double ri3 = ri*ri*ri;
	double ri6 = ri3*ri3;

	force = 48*eps*( pow(sigma,12)*ri6 - pow(sigma,6)/2 )*ri6*ri*ri;
	potential = 4*eps*ri6*( ri6* pow(sigma,12) - pow(sigma,6) );
}

void set_quadres( int &width, int &height, double opt_size )
{
	width = (int)ceil( Lx/opt_size );
	height = (int)ceil( Ly/opt_size );
	cout << "Count of quadres: " << width*height << "\n";
}

int getNumbArea(double x, double y)
{
    int nX = (int) floor(x / (Lx / widthq));
    int nY = (int) floor(y / (Ly / heightq));
    int part = nX*heightq + nY;

    return part;	
}

double optimal_size_of_quadre()
{
	double force, oldforce, potential;
	double r_min = sigma*pow(2.,(1.)/6);
	double size = r_min;
	double step = 0.25*r_min;

	f( size, oldforce, potential );
	size += step;
	f( size, force, potential );

	double prir = (oldforce - force)/abs(oldforce);

	while( prir > 0 )
	{
		size += step;
		oldforce = force;
		f( size, force, potential );
		prir = (oldforce - force)/abs(oldforce);
	}

	while( abs(force) > abs(MISSTAKE*oldforce) )
	{
		size += step;
		f( size, force, potential );
	}

	return size;
}

void computeStartIndexInPart()
{
	startIndexInPart = new int[ widthq*heightq + 1 ];
	int* numberItemsInPart = new int[ widthq*heightq ];

	for( int i = 0; i < widthq*heightq; i++ )
	{
		startIndexInPart[i] = 0;
		numberItemsInPart[i] = 0;
	}
	startIndexInPart[widthq*heightq] = 0;

	for ( int i = 0; i < N; i++ )
	{
		int ind = getNumbArea(x[i],y[i]);
		numberItemsInPart[ ind ]++;
	}

	startIndexInPart[0] = 0;
	for ( int i = 0; i < widthq*heightq; i++ )
	{
		startIndexInPart[i+1] = startIndexInPart[i] + numberItemsInPart[i];
	}
	delete[] numberItemsInPart;
}

void exchangeItem(int _n, int _m)
{
    double tmp;

    tmp = x[_n];
    x[_n] = x[_m];
    x[_m] = tmp;

    tmp = y[_n];
    y[_n] = y[_m];
    y[_m] = tmp;

    tmp = z[_n];
    z[_n] = z[_m];
    z[_m] = tmp;

    tmp = vx[_n];
    vx[_n] = vx[_m];
    vx[_m] = tmp;

    tmp = vy[_n];
    vy[_n] = vy[_m];
    vy[_m] = tmp;

    tmp = vz[_n];
    vz[_n] = vz[_m];
    vz[_m] = tmp;

    tmp = ax[_n];
    ax[_n] = ax[_m];
    ax[_m] = tmp;

    tmp = ay[_n];
    ay[_n] = ay[_m];
    ay[_m] = tmp;

    tmp = az[_n];
    az[_n] = az[_m];
    az[_m] = tmp;
}

void refresh_quadres()
{
	delete[] startIndexInPart;
	computeStartIndexInPart();

	int* curIndexInPart = new int[ widthq*heightq + 1 ];
	for( int i = 0; i < widthq*heightq + 1; i++ )
	{
		curIndexInPart[i] = startIndexInPart[i];
	}

	int j = 0;
	while( j < N )
	{
		int newPart = getNumbArea(x[j], y[j]);
		if( (j < startIndexInPart[newPart+1]) && (j >= startIndexInPart[newPart]) )
			j++;
		else
		{
			int ind = curIndexInPart[newPart]++;
			exchangeItem(j, ind);
		}
	}
	delete[] curIndexInPart;
}

void SortArea()
{
	double* x_tmp = new double[N];
	double* y_tmp = new double[N];
	double* z_tmp = new double[N];
	int* lastIndexInPart = new int[widthq*heightq];

	for (int i = 0; i < N; ++i)
	{
       x_tmp[i] = x[i];
       y_tmp[i] = y[i];
       z_tmp[i] = z[i];
    }

	for( int i = 0; i < widthq*heightq; i++ )
	{
		lastIndexInPart[i] = startIndexInPart[i];
	}

	for( int i = 0; i < N; i++ )
	{
		int ind = getNumbArea(x_tmp[i],y_tmp[i]);
		int newPart = lastIndexInPart[ ind ]++;
		x[newPart] = x_tmp[i];
		y[newPart] = y_tmp[i];
		z[newPart] = z_tmp[i];
	}

	delete[] x_tmp;
	delete[] y_tmp;
	delete[] z_tmp;
	delete[] lastIndexInPart;
}

void init( double _eps, double _sigma, double _dt, ifstream &is )
{
	dt = _dt;
	is >> N ;
	is >> Lx;
	is >> Ly;
	is >> Lz;
	eps = _eps;
	sigma = _sigma;

	double sizeq = optimal_size_of_quadre();
	cout << "\n\nOptimal size of quadre is "<< sizeq << "\n";
	set_quadres(widthq, heightq, sizeq);

	x = new double[N];
	y = new double[N];
	z = new double[N];

	xpr = new double[N];
	ypr = new double[N];
	zpr = new double[N];

	ax = new double[N];
	ay = new double[N];
	az = new double[N];

	vx = new double[N];
	vy = new double[N];
	vz = new double[N];

	for ( int i = 0; i < N; i++ )
	{
		is >> x[i];
		is >> y[i];
		is >> z[i];
		is >> vx[i];
		is >> vy[i];
		is >> vz[i];
	}

	for ( int i = 0; i < N; i++ )
	{
		ax[i] = 0;
		ay[i] = 0;
		az[i] = 0;
	}

	computeStartIndexInPart();
	SortArea();
}

void release()
{
	delete[] x;
	delete[] y;
	delete[] z;

	delete[] xpr;
	delete[] ypr;
	delete[] zpr;

	delete[] ax;
	delete[] ay;
	delete[] az;

	delete[] vx;
	delete[] vy;
	delete[] vz;

	delete[] startIndexInPart;
}

void periodic( double &x, double &y, double &z, double &vz, double &z_prev )
{
	int zCrossingNumber;	
    x = fmod(x, Lx);
    y = fmod(y, Ly);

    if (x < 0)
        x = x + Lx;

    if (y < 0)
        y = y + Ly;

    if( z > Lz ){
        zCrossingNumber = (int) floor(z / Lz);
        if( zCrossingNumber % 2 ){
            z = Lz - fmod(z, Lz);
            vz = - vz;
            z_prev = zCrossingNumber * Lz + (Lz - z_prev);
        }
        else{
            z = fmod(z, Lz);
            vz = vz;
            z_prev = - zCrossingNumber * Lz + z_prev;
        }
    }
    if( z < 0 ){
        zCrossingNumber = (int) floor(abs(z) / Lz) + 1;
        if( zCrossingNumber % 2 ){
            z = - fmod(z, Lz);
            vz = - vz;
            z_prev = - (zCrossingNumber - 1) * Lz - z_prev;
        }
        else{
            z = Lz + fmod(z, Lz);
            vz = vz;
            z_prev = zCrossingNumber * Lz + z_prev;
        }
    }
}

void computing_accs( int i, int j )
{
	double dx = x[i] - x[j];
	double dy = y[i] - y[j];
	double dz = z[i] - z[j];

	if ( fabs(dx) > .5*Lx )
	{
		if ( dx > 0 )
			dx = dx - Lx;
		else
			dx = dx + Lx;
	}

	if ( fabs(dy) > .5*Ly )
	{
		if ( dy > 0 )
			dy = dy - Ly;
		else
			dy = dy + Ly;
	}

	double r = sqrt( dx*dx + dy*dy + dz*dz );

	double force, potential;
	force = potential = 0;

	if( r == 0 )
		r = 1E-15;

	f( r, force, potential );

	ax[i] = ax[i] + force*dx;
	ay[i] = ay[i] + force*dy;
	az[i] = az[i] + force*dz;

	/*ax[j] = ax[j] - force*dx;
	ay[j] = ay[j] - force*dy;
	az[j] = az[j] - force*dz;*/

	P = P + potential;
}

void verlet( int st, int fn )
{
	K = 0;

	for ( int i = st; i < fn; i++ )
	{
		double xnew = 2 * x[i] - xpr[i] + ax[i] * dt * dt;
		double ynew = 2 * y[i] - ypr[i] + ay[i] * dt * dt;
		double znew = 2 * z[i] - zpr[i] + az[i] * dt * dt;

		periodic(xnew,ynew,znew,vz[i], z[i]);

		xpr[i] = xnew;
		ypr[i] = ynew;
		zpr[i] = znew;
	}
}

void periodicArea(int &_x, int &_y)
{
    if (_x < 0)
        _x = _x + widthq;
    if (_x >= widthq)
        _x = _x - widthq;
    if (_y < 0)
        _y = _y + heightq;
    if (_y >= heightq)
        _y = _y - heightq;
}

void accel()
{
	P = 0;
	for( int i = 0; i < N; i++ )
	{
		ax[i] = 0;
		ay[i] = 0;
		az[i] = 0;
	}

	for( int l = 0; l < widthq; l++ )
	{
		for( int k = 0; k < heightq; k++ )
		{
			int currentPart = l*heightq + k;
			for (int i = startIndexInPart[currentPart]; i < startIndexInPart[currentPart+1]; i++)
			{
				for (int n = l - 1; n <= l + 1; n++)
				{
					for (int m = k - 1; m <= k + 1; m++)
					{
						int nX = n; 
						int nY = m;
						periodicArea(nX, nY);
						int neighbourPart = nX*heightq + nY;

						#pragma omp parallel for
						for(int j = startIndexInPart[neighbourPart]; j < startIndexInPart[neighbourPart+1]; j++)
                        {
                            if( (i == j) && (l==n) && (k==m) )
                                continue;

                            computing_accs(i,j);
						}
					}
				}
			}
			verlet( startIndexInPart[currentPart], startIndexInPart[currentPart+1] );
		}
	}

	double *xbuf, *ybuf, *zbuf;

	xbuf = xpr;
	xpr = x;
	x = xbuf;

	ybuf = ypr;
	ypr = y;
	y = ybuf;

	zbuf = zpr;
	zpr = z;
	z = zbuf;

	refresh_quadres();
}

void eiler()
{
	double *xbuf, *ybuf, *zbuf;

	xbuf = x; ybuf = y;
	x = xpr;  y = ypr;
	xpr = xbuf; ypr = ybuf;
	zbuf = z; z = zpr; zpr = zbuf;

#pragma omp parallel for
	for ( int i = 0; i < N; i++ )
	{
		double xnew = xpr[i] + vx[i]*dt;
		double ynew = ypr[i] + vy[i]*dt;
		double znew = zpr[i] + vz[i]*dt;

		periodic(xnew,ynew,znew,vz[i],zpr[i]);

		x[i] = xnew;
		y[i] = ynew;
		z[i] = znew;
	}

	refresh_quadres();
}

void write_molecules( ofstream &os )
{
	for (int i = 0; i < N; i++)
	{
		os << x[i] << " " << y[i] << " " << z[i] << "\n";
	}

}
void write_energies( ofstream &os )
{
	E = K + P;
	os << "\n";
	os << K << " " << P << " " << E;
	os << "\n";
}

void main()
{
	int start = GetTickCount();

	int cProcs = omp_get_max_threads();
	omp_set_num_threads( cProcs );

	printf("Count of threads is: %d\n", cProcs);
	
	ifstream is;
	ofstream os;

	os.open("File.txt");
	is.open("Input.txt");

	double _eps, _sigma; 
	double _dt;

	cout << "dt: ";
	cin >> _dt;

	int determs;

	cout << "determs: ";
	cin >> determs;

	bool autos = true;

	cout << "The type of sigma and epsilon. 1 - Argon, 0 - Inputed value. Your choise:";
	cin >> autos;
	
	if( autos )
	{
		_eps = 165.3E-23;
		_sigma = 3.405;
	}
	else
	{
		cout << "EPSILON = ";
		cin >> _eps;

		cout << "SIGMA = ";
		cin >> _sigma;
	}

	init(_eps, _sigma, _dt, is);

	os << Lx << "\n";
	os << Ly << "\n";
	os << Lz << "\n";
	os << N << "\n";
	os << determs+1 << "\n";

	P = 0;
	K = 0;
	E = 0;

//	write_molecules(os);

	eiler();

//	write_molecules(os);

	for ( int i = 0; i < determs; i++ )
	{
		accel();
	//	write_molecules(os);
	}

	write_molecules(os);
	release();

	os.close();

	int finish = GetTickCount();
	printf("The time of computing is %f\n", finish - start);
}